Pair Correlation Functions and a Free-Energy Functional for the 

Nematic Phase 



Pankaj Mishra, Swarn Lata Singh, Jokhan Ram and Yashwant Singh 
Department of Physics, Banaras Hindu University, Varanasi-221 005, India 

(Dated: February 1, 2008) 

Abstract 

In this paper we have presented the calculation of pair correlation functions in a nematic phase 
for a model of spherical particles with the long-range anisotropic interaction from the mean spheri- 
cal approximation(MSA) and the Percus-Yevick (FY) integral equation theories. The results found 
from the MSA theory have been compared with those found analytically by Holovko and Sokolovska 
(J. Mol. Liq. 82, 161(1999)). A free energy functional which involves both the symmetry con- 
serving and symmetry broken parts of the direct pair correlation function has been used to study 
the properties of the nematic phase. We have also examined the possibility of constructing a free 
energy functional with the direct pair correlation function which includes only the principal order 
parameter of the ordered phase and found that the resulting functional gives results that are in 
good agreement with the original functional. The isotropic-nematic transition has been located 
using the grand thermodynamic potential. The FY theory has been found to give nematic phase 
with pair correlation function harmonic coefficients having all the desired features. In a nematic 
phase the harmonic coefficient of the total pair correlation function /i(xi,X2) connected with the 
correlations of the director transverse fluctuations should develop a long-range tail. This feature 
has been found in both the MSA and FY theories. 

PACS numbers: 71.15.Mb, 64.70.Md, 61.30.Cz 
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I. INTRODUCTION 



The distribution of molecules in a classical system can adequately be described by one 
and two-particle density distributions. The one particle density distribution, p(x) defined 
as 

p(x) = p(r, n) <^(r - r,)S{n - n,) > (1.1) 

i=l 

where Xj indicates both position rj and orientation fij of i^^ molecule, the angular bracket 
represents the ensemble average and S the Dirac function, is constant independent of position 
and orientation for an isotropic fluid but contains most of the structural informations of or- 
dered phases like crystalline solids and liquid crystals. The two-particle density distribution 
p(xi,X2) which gives probability of flnding simultaneously a molecule in volume element 
dridfli centered at (ri,r2i) and a second molecule in volume element dr2dQ2 centered at 
(r2, ^2) is defined as 

p(xi, X2) = p(ri, ill, r2, n^) =< J2 '^(ri - r^)6{n, - n,)6{r2 - r^)6{n2 - ^j) > (1.2) 

i¥=3 

The pair correlation function g'(xi,X2) is related to p(xi,X2) by the relation 

P(X1,X2) , . 

ff(Xi,X2) ^ (1.3) 
p(Xi)p(x2) 

Since in an isotropic fluid p(xi) = p{^2) — Pi — where < > is the average number 
of molecules in volume V, 

p^g{r,n,,Q2)=p{r,n,,n2) (1.4) 

where r = (r2 — ri). In the isotropic fluid gr(xi,X2) depends only on inter particle distance 
I r2 — I'll — r, orientation of molecules with respect to each other and on the direction of 
vector r (f = Ms a unit vector along r). These simphflcations are due to homogeneity which 
implies continuous translational symmetry and isotropy which implies continuous rotational 
symmetry. Such simplifications do not generally occur in ordered phases. 

The pair correlation functions as a function of intermolecular separations and orientations 
at a given temperature and pressure can be found either by computer simulation [1-5] or by 
simultaneous solution of an integral equation, the Ornstein-Zernike (OZ) equation, 

/l(Xi,X2) = c(Xi,X2) + j c(Xi,X3)p(x3)/l(x3,X2)(iX3 (1.5) 
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where = dr^dQ,^ and /i(xi,X2)(= 5'(xi,X2) — 1) and c(xi,X2) arc respectively, the total 
and direct pair correlation functions, and an algebraic closure relation which relates the 
correlation functions to the pair potential. Well known approximations to the closure relation 
are the hypernetted-chain relation, the Percus-Yevick (PY) relation and the mean spherical 
approximation (MSA) [6]. These integral equation theories have been quite successful in 
describing the structure and thermodynamic properties of isotropic fluids [7-11]. However, 
their application to ordered phases which can be regarded as inhomogeneous, have so far 
been very limited [12-15], though no feature of the theory inherently prevents them from 
being used to describe the structure of ordered phases. One of the problems that arises in the 
case of ordered phases is the appearance of p(r, Q) in the OZ equation (see Eq.(1.5)). This 
implies that in contrast to the isotropic case where we needed only two relations, namely the 
OZ equation and a closure relation, an additional relation corresponding to single particle 
distribution connecting to pair correlation function is needed to solve the ensuing equations 
self consistently. 

In this paper we take nematic in which molecules are aligned on the average along a par- 
ticular but arbitrary direction while the translational degrees of freedom remain disordered 
as in an isotropic phase, as an example of an ordered phase. At the isotropic-nematic transi- 
tion the isotropy of the space is spontaneously broken and as a consequence, the correlations 
in the distribution of molecules lose their rotational invariance. The change from isotropic 
fluid to nematic state in the absence of external field involves collective fiuctuations, which 
develops orientational wave excitations known as Goldstonc modes[16]. This leads to the 
divergence of the corresponding harmonics of the total pair correlation function /i(xi,X2) in 
the limit of zero wave vector. By computer simulation of a system of ellipsoids Phoung and 
Schmid [13] have evaluated the effect of breaking of rotational symmetry on pair correlation 
functions and showed that in a nematic phase there are two qualitatively different contribu- 
tions; one that preserves rotational invariance and the other that breaks it and vanishes in 
the isotropic phase. 

Holovko and Sokolovska [14] have used the MSA closure relation and the Lovett equation 
[17] (see Eq.(3.17))which relates one particle density to pair correlation function to solve 
analytically the OZ equation for a model of spherical particles with the long range anisotropic 
interaction (see Eq.(2.1)) in a nematic phase. However, when Phoung and Schmid [13] used 
the PY closure and the Lovett equation and solved the OZ equation numerically for a 
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system of soft ellipsoids, nematic phase was not found and for this the PY closure was 
blamed. Zhong and Petschek [18] have analyzed the diagrammatic expansion of the direct 
correlation function and concluded that the PY closure can not reproduce the Goldstone 
modes in the general case of spontaneous partial ordering. 

Recently we [19] used the PY closure and solved numerically the OZ equation for a sys- 
tem of elongated rigid molecules interacting via the Gay-Berne potential [20] and showed 
that the PY closure gives nematic phase with the pair correlation function harmonic coeffi- 
cients having features similar to those found by computer simulation [13] and by analytical 
solution [14]. Instead of using a closure relation for p(r, ft) we expressed it in terms of order 
parameters and solved the resulting equation for values of order parameters ranging from 
zero to some maximum value. Non-zero values of order parameters break the symmetry of 
isotropic phase and the degree of symmetry breaking is given by the values of the order 
parameters. Using these correlation functions we constructed a free energy functional and 
used it to determine the value of order parameters in the nematic phase by minimizing it. 
Once the values of order parameters are known the pair correlation functions in the nematic 
phase are obtained from the known results. 

In this paper we extend our method to calculate the pair correlation functions in nematic 
phase using the MSA and PY closure relations for a system the molecules of which interact 
via a pair potential considered in ref. [14]. This allows us to compare our results for the MSA 
with those found analytically and therefore to test the accuracy of our method. The PY 
relation is shown to give nematic phase with all the expected features. The paper is organized 
as follows: In Sec.II we describe the MSA and PY integral equation theories and give a brief 
account of computational procedure. In Sec. Ill we construct a free energy functional of an 
inhomogeneous system that contains both symmetry conserved and symmetry broken parts 
of the direct pair correlation function. The isotropic-nematic transition point and freezing 
parameters are calculated in Sec IV. The paper ends with discussions given in Sec. V. 

II. CORRELATION FUNCTIONS 

The pair potential used by Holovko and Sokolovska [14] in their analytical solution of the 
OZ and Lovett equation with MSA closure has the form 

^;(xi, X2) = Vhs{r) + vo{r) + V2{r, l^i, Q2) (2.1) 
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where Vhs is the hard sphere potential 



Vhs{r) =00 r < a (2.2) 
= r> a (2.3) 



The long-range attraction has isotropic part 



2exp(-2;or) 



^^o(r) = -aoizoay " ^ (2.4) 

r/a 



and the anisotropic part 



2exp(-Z2r) 



f2(r, fii, 1^2) = -a2(-22cr) -/ P2(cosr2i2) (2.5) 

r/(7 

where P2(cos ^,12) is the second order Legendre polynomial of relative molecular orientations. 
This model potential is independent of orientation of the intermolecular separation vector r. 
This fact limits the number of harmonic coefficients that appear in the spherical harmonic 
expansion of pair correlation functions. 

We choose a coordinate frame with it's z-axis in the direction of the director ^(director 
frame). The director n is a unit vector along the direction of ahgnment of molecules. All 
orientation dependent functions are expanded in spherical harmonics y;OT(Q)[6]. This yields 
(for uniaxial nematic phase of axially symmetric molecules) [21] 



p 



4^ Z(even) 



E /^^'o(^^) (2.6) 



where /; = y (2/ + l)Pi and /o = 1. Pi for / 7^ are order parameters; their values are 
zero in the isotropic phase and nonzero in the nematic phase. For two particle functions one 
has [22] 

^(r, 1^1,^2)= E i^hi2irmm2m{r)Yi,mA^i)Yi,^,{n2)Yil{r) (2.7) 

Iil2lmim2m 

where ip stands for /i or c or v. In uniaxial nematic phases, only real coefficients with 
777.1 + 7712 — m = and even /1 + /2 + / enter in the expansion. Since the molecules in the model 
system under consideration have axial symmetry, every single I is even as well. Because, in 
isotropic phase h and c preserve the rotational symmetry, for them 

'iphi2imim2m{r) = i>hi2i{r)Cg{hl2lmim2m) (2.8) 



where Cg is the Clebsch-Gordan(CG) coefScient. The absence of the CG coefficients 
in Eq.(2.7) when ip represents pair correlation functions of nematic phase, removes the 
restriction \li — < I < h + I2 on the values of the index I. As a consequence, coefficients 
such as ■0200000 (j") and ■0o2oooo('") are nonzero in nematic whereas they do not survive in the 
isotropic case. The emergence of these harmonic coefficients are due to symmetry breaking. 

To solve the OZ equation it is advisable to use the Fourier representation. The expansion 
coefficients 4'hi2iraim2m{f) are related to their counterparts in Fourier space by the Hankel 
transform 

'4^hi2imim2m{k) = 47ri' / drr"^ ji{kr)^i^i^im^^^m{r) , (2.9) 

Jo 

47r( i)' 

1phl2lmim2m{r) = . .3 / dkk'^ ji{kr)'^i^i^imim2m{k) (2.10) 

where ji{kr) is the spherical Bessel function. 

Using Eqs(2.9)-(2.10) and the spherical harmonic expansion for the correlation functions 
(Eqs.(2.6) and (2.7)), OZ equation reduces in the k-space to the form 

'Hil2lmiTn2m ^I\l2lm\m2m 

{k) 

^hl2l'mim2m 

{k) 



mim„m 

\k)h, 

0I2I m„m2m 



/47r Y 

T'r% r", (2.11) 

rri rn n m m m V / 



where p* = pa^ , the symbol indicates the collection of nine indices ,m'^, 
m'^,m',m" and notation 



^'r^!L2m = / dnYCMyhmMyi2m2{^) 



\\ ^^\t(2/+lf Cg{hl2mO)Cg{hymim2m) (2.12) 

Since the pair potential of Eq(2.1) is independent of orientation of the intermolecular 
separation vector r, the harmonic coefficients that survive in the expansion of pair correlation 
functions have I — m — and mi = —1712- This allows notational simplification from six 
indices to three. We therefore rewrite Eq. (2.11) as 

lhl20m-mo{k) = Ihhmik) 



£ E %i',Jk)h,,^i2M Y: V2rTTrS (2.13) 

'3*3 



where m = —m. 

Since there is no summation over index m on the right hand side of Eq.(2.13), the OZ 
equation for harmonics with different values of m decouple. The equation corresponding to 
the isotropic case is found by putting Po — 1 and Pl^q = in Eq(2.13). 
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FIG. 1: Comparison of some of the harmonics of the total pair correlation function in the Fourier 
space for the nematic phase {(3a2 = l,/3ao = 0.1, zoa = z^o = 1, rj = 0.315, P2 = 0.63, P4 = 0.27) 
obtained by the analytical solution[14] (dashed line) with the results found using numerical method 
(full line) for the MSA. The two curves are indistinguishable at the scale of the figure. 



A. MSA Closure 

The MSA relation is written as 

/i(xi,X2) = -1, |r2-ri|<(j (2.14) 
c(xi,X2) = -/3T;(r,Qi,Q2), |r2-ri|>(7 (2.15) 

where v is given by Eqs.(2.4) and (2.5). 
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FIG. 2: Comparison of the harmonic coefficients C22i{t*) and /i22i(''*) obtained by the analytical 
solution[14] (dashed line) with the results found using analytical method (full line) for the MSA. 
Parameters are same as in Figure 1. A 1/r* tail in harmonic coefficient ^221 (^*) is seen. 



Condition (2.14) is exact for potential model of Eq(2.1) since g{r, fli, 0,2) = for r < a. 
However, it is only for large r that c(r, 0,2) is asymptotic to (3v{r, fli, fl2) but in MSA it 
is assumed that c = —f^v for all r(r > a). 

Using Eq.(2.7) with the condition Z = m = to expand the pair correlation functions 
and Eq(2.14) we get for r < cr 



cooo(r) = -7ooo(r) - (^nf^ (2.16) 

and 

Chhmir) = -7hi2m(r) (2.17) 
For r >(7 from Eqs.(2.15), (2.4) and (2.5) we get 

Cnm[r)-[ ij + i r/a ^ ^ 

where i — 0,2 and m — 0,1, 2. 




k/2n 

FIG. 3: Comparison of the structure factor curves for the nematic phase (/3a2 = 1, /3ao = 0.1, zqct = 
2:20" = 1, 77 = 0.315, P2 = 0.63, P4 = 0.27). The dashed hue is obtained with the analytical 
solution [14] while full line is obtained by our numerical method for the MSA. The two curves 
overlap at all values of k. 

B. The PY Closure 

The PY relation is written as 



c(Xi,X2) = /(xi,X2)[c/(Xi,X2) -c(Xi,X2)] 



(2.19) 



where /(xi, X2) = exp[— /3^;(xi, X2)] — 1 is the Mayer function and /3 = (/csT)"^; ks being the 
Boltzmann constant and T, temperature. Expansion in spherical harmonics with constraint 
I — m = leads to 



1 



,/ / 1,1'," II 



/ // 
m m m 



1 / // 

m m m 



(2.20) 



where ff^i'^fir) is the harmonic coefficient of the Mayor function /(r, Jli, 0,2)- 

We solved the OZ equation for both the MSA and the PY closures for given values of 

order parameters P2 and P4. In order to solve these equations numerically we followed the 
iterative method described in ref.[13]. However, as coefficients /ijjjjml^) decay slowly 
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FIG. 4: Plot of the harmonic coefficients of the direct pair correlation function C22o{r*), C22i{r*) 
and C2oo{r*) obtained by solving PY integral equation theory (/3a2 = l,/3ao = 0.1, zga = Z20 = 1, 
rj = 0.30). While the contribution of symmetry breaking part in C22i(r*) is small (shown by dot- 
dashed line), in C22i(r*) it is comparable inside the core. The harmonic coefficient C22o(^*) arises 
due to symmetry breaking only. 

and extend to relatively large values of r*(= ^) in nematics we have extended the range of 
r*(i.e. r* = 60) to ensure proper convergence. The other point which needed special care is 
related to the pronounced long range tail which occurs in coefficients h\^i^^{r*) with m = ±1 
(see Fig 6). Before performing the Hankel transform in each iteration we fit the data points 
of these harmonic coefficients beyond r* > rQ{— 20) to a power law a + ^, shift them by a 
and then extrapolate them to infinity [19]. This removes the finite size effect on the tail. 

The potential parameters taken in our calculations are zoa = Z2(j = 1, /3ao = 0.1 and 
(3a2 = 1. For the PY wc have also considered the case of f3a2 = 0.5. 

In Fig.l wc compare the results of the Fourier transform of some of the harmonic co- 
efficients hi^i^m{r*) obtained with the analytical solution of the model potential [14] with 
the results found using numerical method stated above for the MSA. Both results are for 
r](— ^) = 0.315, P2 — 0.63 and P4 — 0.27. The minimization of free energy functional 
(see Sec III) which contains both the symmetry breaking and symmetry conserving parts 
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of the direct pair correlation function gives these values of order parameters P2 and P4. 
These values of order parameters are also found from the analytical result of ref[14]. Both 
curves shown in the figure overlap indicating an excellent agreement between the two re- 
sults. In Fig. 2 we compare the harmonic coefficients C22i(r*) and h22i{r*). These harmonic 
coefficients are of fundamental importance as they appear in nematic elastic constants. The 
decay of h22i{r*) as ^ at large distance is clearly seen. In Fig. 3 we compare our results of 
the structure factor defined as 



Again we find excellent agreement between analytical and numerical results including 
small peak at k=0 which is attributed to the appearance of additional effective attraction 
due to parallel alignment of molecules [14]. 

In Figs 4-6 we give results found from using the PY closure for /?a2 = l,r] — 0.30, P2 = 
0.69 and P4 — 0.32 in the director space. These values of order parameters have been 
found from the minimization of the free energy functional (see Sec III) . While the harmonic 
coefficients C22o(''*) and C22i(r*) shown in Fig 4 survive both in the isotropic P2 — P4 — 
and in the nematic phase P2 7^ 0,P4 7^ 0, the harmonic coefficient C2oo{r*) survive only 
in the nematic phase and vanishes in the isotropic phase. The contribution arising due 
to symmetry breaking to the harmonic coefficients C22o(?^*) and C22i(r*) shown in Fig 4 by 
dot-dashed line are found to be very small compared to the symmetry conserving part. Few 
selected harmonic coefficients of h are shown in Figs 5 and 6 in the director space. In Fig 

5 we plot the harmonic coefficients /i2oo(^'*) and h22o{r*). While h2m{^*) survive only in 
the nematic phase, /i22o(^*) survive both in the isotropic and in the nematic. In the case 
of /i22o(^*) we also plot the contributions arising due to symmetry breaking and symmetry 
conserving and note that the contribution arising due to symmetry breaking is small. In Fig 

6 we plot harmonic coefficients /i22i('"*) and show its ^ dependence in the inset. 

In Figs 7 and 8 we plot few selected harmonic coefficients of c and h in director space for 
/3a2 = 0.5, 77 = 0.48, P2 — 0.54 and P4 = 0.12. As will be shown later that at the isotropic- 
nematic transition the packing fraction 77 of the nematic phase is 0.458. for I3a2 = 0.5 while 
at j3a2 = 1 it is 0.244. The comparison of these harmonic coefficients show that orientational 
ordering has more pronounced effect on these harmonic coefficients when (5a2 — l compared 




1 + 



(47r)3/2 



[hmo{k) + 2^/5P2h2oo{k) + hP^h22o{k)]. 



(2.21) 
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to that of /3a2 = 0.5. This could be easily understood from the fact that the orientational 
ordering arises solely due to the long range anisotropic part of the interaction and j3a2 
measures its strength. We note that the PY closure gives harmonic coefficients of both 
symmetry breaking and symmetry conserving parts of pair correlation functions which have 
features similar to those found from the MSA solution as well as from computer simulation 
[13]. 
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FIG. 5: Coefficients ^200 (j'*) and /i22o(''*) obtained from the PY theory with the parameters same 
as in Figure 4. The symmetry breaking contribution shown by dot-dashed fine is very small for 
h22o{T*)- The harmonic coefficient /i20o('^*) arises due to symmetry breaking only. 

In Fig 9 we plot the structure factor found from the PY theory. In this case its expression 
is found to be 

S{k) = 1 + E ^{2h + l){2l, + l)Pi,PiM.Uk) (2.22) 

The curve drawn in full hue corresponds to 77 = 0.30, P2 — 0.69, P4 = 0.32, /3a2 = 1 whereas 
the one drawn in dashed line corresponds to 77 = 0.48, P2 = 0.54, P4 = 0.12 at /3a2 = 0.5. 
We note that 77 = 0.48 is close to the freezing transition where the system goes into the 
crystalline phase and therefore the peaks in S{k) are more pronounced compared to the 
curve corresponding to 77 = 0.30. According to Hansen- Ver let [25] criterion fluid becomes 

12 




-0.24 - 




0.15 



I/r 



-0.30 







3 6 9 



15 



r 



FIG. 6: Harmonic coefficient /i22i('"*) in the director frame. Details are same as in Figure 4. Inset 
shows the plot of h22i{r*) with respect to l/r*; the dashes line shows the extrapolated part. The 
origin of tail is due to orientational symmetry breaking. 

unstable when the height of main peak in S{k) becomes equal to 2.9 ± 0.1. The curves 
corresponding to rj = 0.30, /?a2 = 1 shows a small peak at A; = as was found in the 
case of the MSA theory. However, this peak is not seen in the curve corresponding to 
rj — 0.48, — 0.5. For this case as indicated by the values of order parameters, the 
orientational ordering is weak compared to that of 77 = 0.30, /3a2 = 1 and therefore the 
effective attraction which arises due to orientational ordering is negligible. 

III. FREE ENERGY FUNCTIONAL 

The reduced free energy A[p] of an inhomogeneous system is a functional of p(x) and is 
written as [21] 



A[p] ^ Aid[p] + Aea^ip] 



(3.1) 



The ideal gas part is exactly known and is given as 




(3.2) 
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FIG. 7: Plot of the symmetry breaking harmonic coefficients (c2oo('"*)) ^20o('"*)) and symmetry 
conserving harmonic coefficients (C220 (''*)) ^220 (''*)) found by PY integral equation theory with the 
parameters /?a2 = 0.5, /?ao = 0.1, zoa = z^o = 1, 77 = 0.48. 



where A is cube of the thermal wavelength associated with a molecule. The excess part aris- 
ing due to intermolecular interactions is related to the direct pair correlation function(DPCF) 
as 

^^;(|^ = -c<")(x.,x.;po)-c«»>(x.,x.;lpl) (3.3) 

where superscripts (0) and in) represent respectively the symmetry conserving and symme- 
try breaking parts of the DPCF. In other words, c*^"-* is found by putting order parameters in 
Eqs(2.13) equal to zero whereas c*^") are the contributions which arise when order parameters 
are nonzero. 

is found by functional integration of Eq.(3.3). In this integration the system is taken 
from some initial density to the final density p(x) along a path in the density space; the 
result is independent of the path of integration [26]. For the symmetry conserving part c*^°) 
the integration in density space is done taking isotropic fluid of density pi (the density of 
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FIG. 8: Harmonic coefficients C22i(r*) and h22i{r*) obtained by PY theory with the parameters 
same as in Figure 7. 



coexisting fluid) as reference. This leads to 

AfJ[p] = AM - ^ / ^^i/ t^X2Ap(xi)Ap(x2) 

xc(°)(xi,X2) (3.4) 

where 

c('^)(xi,X2) ^2jd\\j dA'c(°){xi,X2;p^ + AA'(po-pO} 

and 

Ap(x) = p(x)-p^ (3.5) 

Aexi.Pi) is the excess reduced free energy of isotropic fluid of density pi and po is the average 
density of the ordered phase. 

In order to integrate over c*^"^[p], we characterize the density space by two parameters A 
and ^ which vary from to 1[19]. The parameter A raises density from to po as it varies 
from to 1 whereas parameter ^ raises the order parameter from to as it varies from 
to 1. If we have n order parameters to describe the ordered phase we can think of a 
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n— dimensional order parameter space; a point in tliis space defines tiie values of tfie n order 
parameters. Tlie integration over c^'^\p) in Eq(3.3) can be done along a straight line path 
that connects origin to a point corresponding to the final values of all n order parameters. 
This path is characterized by the variable ^. This gives 

4?[P] = -ll d^i J c?X2p(xi)p(x2)S(^)(xi,X2) (3.6) 

where 

c(")(xi,X2) = 4 d^^ C C d\\ dX X 

Jo Jo Jo Jo 

n 

c(")(xi,X2,AA'po;e^\ (3-7) 

\ 1=1 

While integrating over A the order parameters Pi are kept fixed and while integrating 
over ^ the density is kept fixed. The result does not depend on the order of integration. The 
free energy functional of an ordered phase is the sum of Aid, and ^4^"^ given respectively 
by Eqs(3.2), (3.4) and (3.6). Note that the Ramakrishnan and Youssouff [23] free energy 
functional is the sum of only A-i^ and A^J and contains an additional approximation in which 
c*^°^(xi,X2) in (3.4) is replaced by c(xi,X2;pi). 

The minimization of AA — A[p] — A{pq) where A{pq) is the free energy of an isotropic 
phase of density po leads to 

ln/(f]) = C + J dx2Ap(x2)c(°)(xi,X2,po) + / dx2p(x2)ci")(xi,X2) (3.8) 



where /(f2) = and 



) (3.9) 



fl rl rl " 

c^(xi,X2) = 2 / dXX / dX' / dec(")(xi,X2;AA'po;e, E^' 
Jo Jo Jo 

c(°)(xi,X2) = ydAc(°){xi,X2;pi + A(po-pO} 

(3.10) 

The constant C is found from the normalization condition 

J f{fl)dfl = 1 (3.11) 

In order to evaluate c*^"'^(xi, X2) and Ci"^(xi,X2) from Eqs.(3.7) and (3.9) we need symmetry 
breaking part of DPCF from density zero to po and order parameters form zero to Pi at 
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sufficiently small intervals. Tlie computational time needed to evaluate these correlation 
functions depends on the number of order parameters one takes in the calculation. A nematic 
is adequately described by two order parameters, P2 and P4. However, ordered phases such 
as smectic and crystalline solids may need several order parameters. It is therefore advisable 
to approximate the values of A^J^ [p\ with as small number of order parameters as possible. 
Here we show that for nematic it is a good approximation to consider only P2 in calculating 
from Eq.(3.6). 

In case of the MSA the free energy functional reduces to 

4p(x)]-4Po] ^Md^ -lnZ+[P|a^°i(0) + P.4S(0) + P^^^^^^^^ 

-\pA{^) - ^[aS(o) + P2(cti(o) + aS(o)) 

+P|gS(0)]. (3.12) 

where 

Z = i f_^dcos{e) exp[(P2gS(0) + ci2(0) + P2c;2(0))P2(cos 9)] (3.13) 

and 

^S.o(O) = ^V(2^i + 1)(2^2 + 1) [<^ijrydr. 

Note that in this case P4 does not appear exphcitly. It appears only through c^") and c^f^ . 
We calculated c*^"-' and Ci'' a,i r] = 0.315 using the values of the DPCF obtained for r] from 
zero to 0.315 at the interval of 0.02, P2 from zero to 0.70 and P4 from zero to 0.35 at the 
interval of 0.05. Substituting the values of c*^°^ (which correspond to P2 = P4 = 0) , c*^"^ and 
Ci"^ in Eq. (3.12) we minimized the free energy with respect to P2. The order parameter P4 
is found from the relation 

^4 = ^ /cicos(^)exp[(P25S(0) + 512(0) + P2c''iS(0))P2(cos^)]P4(cos^), (3.14) 
The values found are P2 = 0.63, P4 = 0.27. The values of structural parameters defined as 

Q,^,o(0) = -^7(2/i + l)(2/2 + l) r ci,i,o{rydr (3.15) 
^/A■K * Jo 

and 

roo 

Qii2i(0) = Po / ci^i2i{r)r^dr (3.16) 
Jo 

are 6220(0) = 5.29 and C22i(0) -3.60. 
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We next calculate c*^"-* and Ci in same way except taking P4 = 0. When these values of 
c(") and were used we found P2 = 0.65, P4 = 0.28,C22o(0) = 5.46 and C22i(0) = -3.62. 
The two set of values compare well and indicate that using only principal order parameter 
in calculating A'^^ is a good approximation. 

The Ward identity which must be satisfied in a nematic phase relates the single particle 
distribution to an integral of direct pair correlation function. When this identity is expressed 
in a functional differential form it reduces to the Lovett equation [17] 

VnMp{xi) = J c{r,ni,n2)VnA^2)drdn2 (3.17) 

Expanding it in spherical harmonics we get 

1 = -pE\/^^|^(2^' + l)^,'C,(Zir2000)C,(Zir2101) J drr\,2i(r) (3.18) 
III' 

When the two sets of parameters reported above, one in which both P2 and P4 appeared in 
calculating the pair correlation functions while in other only P2 appeared, are substituted in 
this equation we find that it is satisfied with accuracy better than 10~^. From these results 
we conclude that it is sufficient to evaluate c*^"^ and c^^^ with principal order parameter only. 
All results given below correspond to this approximation. 
For the PY the free energy functional is found to be 

A[pi.)]-A[p.] ^AA^ ^ -lnZ+[P|4o(0) + M2(0) + P|a.S(0) + P2P4a.2 

+P!cliO) + P4£ii:o(0) + P2P4c;S(0) + Pf gi2(0)] 

-lipfc^^) + 2P2P44i(o) + prcZm 

—[cSoio) + ^2(^2(0) + aS(o)) + p|sS(o) 

+2P4c;2(0) + 2P2P4a2(0) + P4'g2(0)]. (3.19) 

with 

Z ^ I l\cos{9) exp[{P25S(0) + ^iSW + ^251^(0) + P4CiS(0)}/'2(cos^) + 

{i"4cl4o(0) + di^oi^) + P2Ci2(0) + P4Cu2(0)}P4(cose)] (3.20) 

In these equations, unlike the MSA (see Eq.(3.12)), both P2 and P4 appear. The values 
of c|";^2m(0) ^U^hmi^) have been calculated at r] — 0.30 using the values of the harmonic 
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coefficients of c*^") obtained for r] from zero to 0.30 at the interval of 0.02 and the P2 from 
zero to 0.75 at the interval of 0.05. We substituted these values of c^^^ and c^"^ and in 
Eq(3.19) and minimized the resulting expression with respect to P2 and P4. The values of 
these order parameters and the values of the structural parameters at 77 = 0.30 and /3a2 = 1 
are found to be: 

P2 = 0.69, P4 = 0.32, C22o(0) = 5.64, 6221(0) = -3.84, 6421(0) = 0.045. 
When these values were used in Eq.(3.15) the Ward identity was found to be adequately 
satisfied. 

When we chose /9a2 = 0.50 then as shown in the following section the nematic packing 
fraction at the nematic-isotropic transition is found to be 77 = 0.458. We therefore calculated 
the free energy at 77 = 0.48 which is in the nematic region. In this case we found 

P2 = 0.54, P4 = 0.12, C22o(0) = 4.58, £221(0) = -3.11 ,£421(0) = 0.048. 

These values also satisfy the Ward identity. 

IV. ISOTROPIC-NEMATIC TRANSITION 

The grand thermodynamic potential, defined as 



where is the chemical potential, is generally preferred to locate the freezing transition as it 
ensures that the pressure and chemical potential of the two phases at the transition remains 
equal. Using Eqs.(3.1)-(3.6) and (4.1) wc get 




(4.1) 





Pf 

dxi / (ix2Ap(xi)Ap(x2)c('')(xi,X2) 




(4.2) 



where 



Ap{x) 



P{x) - Pi 

^[(1 + Ap*){l + 5P2P2{cose) + 9P4P4(cos^)} - 1] 



(4.3) 



and 



Ap* 



Po- Pi 



(4.4) 



Pi 
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FIG. 9: Plots of structure factor obtained by using PY theory for (3a2 = 1 (r? = 0.30, P2 = 
0.69, P4 = 0.32) (full line) and /3a2 = 0.5 (r/ = 0.48, P2 = 0.54, P4 = 0.12) (dashed line). Other 
potential parameters are same as in figure 7. A small peak at A; = exists for (3a2 = 1 but not for 
/3a2=0.5. 



The minimization of Al^ with respect to p(x) leads to the following relations for order 
parameters 

J-i Pi{cos 6) cxp{Ii)dcos6 



Pi = 



^\ exp(Ii)dcos6 



where 



and 



1 + Ap* = - / exp(72)rfcos6l 

/i = y dx2Ap(x2)c^°^(xi,X2;po) + j (ix2Ap(x2)cS''^(xi,X2;p) 
h = j (ix2Ap(x2)c(°)(xi,X2;po) + j cix2Ap(x2)4"Hxi, X2; p) 



(4.5) 
(4.6) 

(4.7) 
(4.8) 
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FIG. 10: P* — p* isotherms for the MSA and PY closures for = 1. Other potential parameters 
are same as in Figure 1. The plateau corresponds to change in densiy at the isotropic- nematic 
transition. 

The c1 is defined by Eq(3.9) and 



Jo Jo Jo 



\J:p") (4-9) 

JO JO 

In the isotropic phase the order parameters become zero. Eqs(4.5)-(4.6) are solved self- 
consistently using the values of c^°\ c^"^ c^^\ Ci"^ and c^^^ evaluated in the previous sections. 
By substituting these solutions in the expression of AW we locate the transition. At a given 
temperature and density a phase with lowest grand potential is taken as the stable phase. 
Phase coexistence occurs at the value of pi that makes —AW/N = for the nematic and 
liquid phases. The results arc given in Table 1 for both the MSA and PY theories. 

The pressure can be found using the compressibility equation. In the case of isotropic 
phase 

— = 1-- rdpcooo{0:P) (4.10) 
p p Jo 



p p 

For the nematic phase the relation is found to be 
i3P 1 fp dp 



P ^0 1 + (4372 / dvT^ V(2^i + 1)(2^2 + l)Pi,Pi,hi,Ur, p') 
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(4.11) 




0.65 0.70 0.75 0.80 0.85 0.90 0.95 

* 

P 

FIG. 11: P* — p* isotherm obtained for the PY closure with potential parameter ^a^ = 0.5. Other 
potential parameters are same as in figure 1. 

TABLE I: Isotropic-ncmatic transition parameters of the model potential with /3a2 = 1 and f]a2 = 
0.5 keeping other parameters fixed at zqct = Z2cr = l,/3ao = 0.1. The pressure is given as P* = 
/3P/p. 



P0 2 


Closure 


PI 


Ap* 


Pi 


Pi 


P* 


1.0 


MSA 


0.490 


0.120 


0.540 


0.150 


1.480 


1.0 


PY 


0.411 


0.134 


0.450 


0.130 


1.139 


0.5 


PY 


0.864 


0.013 


0.440 


0.120 


7.700 



In Fig 10 we compare the pressure found from the MSA and PY theories for f3a2 = 1. In Fig 
11 the pressure found from the PY theory is given for /3a2 = 0.5. The plateau corresponds 
to the change in density at the transition. The value of pressure found at the transition is 
given in Table 1. 
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V. DISCUSSIONS 



We have presented the calculation of the pair correlation functions ^(xi, X2) and c(xi, X2) 
in a nematic phase for a model of spherical particles with the long range anisotropic inter- 
action from the MSA and the PY integral equation theories. We chose this model system 
because for this the OZ equation with the MSA closure has been solved analytically [14]. 
The inhomogeneous OZ equation involves the single particle density distribution p(x) which 
we have expressed in terms of order parameters. Non-zero values of order parameters break 
the rotational symmetry. The value of order parameters determine the degree of symmetry 
breaking. For determining the value of order parameters at the isotropic-nematic transition 
we used the equation which we found by minimizing the grand thermodynamic potential. 
The transition from the isotropic to nematic in the density-temperature plane is found by 
solving simultaneously the equations for the grand thermodynamic potential and the order 
parameters. This solution gave the value of density, temperature and order parameters at 
the transition which we have listed in Table 1. The value of order parameters in the nematic 
region has been found by minimization of the reduced Helmholtz free energy functional in 
terms of order parameters. 

The free energy functional given here includes both the symmetry conserving and sym- 
metry broken parts of the DPCF and therefore correctly describes the ordered phase. In the 
free energy functional of Ramakrishnan and Yousouff[23] the DPCF of the ordered phase is 
replaced by that of the coexisting isotropic fluid. This amounts to neglecting the symmetry 
breaking part of the DPCF. In the weighted-density approximation of Curtin and Ashcroft 
and various versions of it [24] the free energy functional is constructed in such a way that 
the free energy density of an inhomogeneous system at a given point is replaced by that 
of a homogeneous system but taken at an auxiliary density which depends parametrically 
on the chosen point. This approach also neglects the new features that emerge in the pair 
correlation functions due to symmetry breaking. 

One of the important features of the total pair correlation function of a nematic phase is 
the appearance of 1/r* tail in harmonic coefficient /iziZ20mii-io(or in om notation hi-i^i^i{r*)) of 
/i(xi,X2). This has been seen in the computer simulation [13] and in the analytical solution 
of the MSA theory by Holovko and Sokolovska[14]. We found this feature in both the MSA 
and PY theories. The long-range tail behaviour of hiii^i{r*) is attributed to the director 
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transverse fluctuations which give rise to the Goldstone modes. This can be seen by taking 
the tensor order parameter Qq,^ = J2i^i ^{^ia^iiJ — '^Said) where a, P = x,y, z and Cj^ is the a 
component of the molecular axis vector ei of each molecule and Sa/s the Kronecker symbol and 
calculating (assuming that the director is along z and the y— axis is perpendicular to wave 
vector k) the correlation {Qxzik)Qxzi—k)) ■ The result involves coefficients ^[";2«mim2m(^) 
with \mi\, \m2\—l. These coefficients which are the Fourier transform of hl^^limim2mi''^*) 
behave as for k ^ 0. 

The calculation of A^^-'fp] involves integration over c*^"-* (xi , X2 ; [p]) in the density space. 
This space is characterized by two variables A and ^ which vary from to 1 and raise the 
density from to po (the average number density of the ordered phase) and order parameters 
from to their final values. One has therefore to evaluate c^")(xi, X2; [p]) from zero value of 
number density and order parameters to their final values at small intervals. This may need 
large computational investment. We have therefore examined the possibility of calculating 
by considering the DPCF which involves only the principal order parameter i.e. P2. 
The resulting free energy functional has been found to give (see Sec III) results which are 
close to the exact results. 

It is important to note that the density-functional approach allows one to include more 
order parameters in the theory even though they are not included in calculating A^") [p] . This 
is done through the parametrization of /)(x)[21]. The results given in Table 1 and in Sec III 
for the PY theory correspond to this approximation. The harmonic coefficients plotted in 
Figs 1-9 have been calculated using both P2 and P4. 

The harmonic coefficient hi^i2i{r*) has been found to be sensitive to values of P2 and P4. 
While P2 creates the 1/r* tail (or equivalently makes its Fourier transform to diverge at 
k=0) P4 suppresses it. This can be seen from the OZ expression 

^"^^^^ = l+(4^[l + 0.7T4p!-1.714P4]a22i(fc) ^'-'^ 

When at A; — > the second term in the denominator approaches to -1 the divergence 
occurs. Since C22i{k — > 0) is negative the term involving P2 help while the term involving 
P4 opposes the divergence. 

The theory developed here can be extended to other ordered phases. Since the symmetry 
breaking part of pair correlation functions have features of the ordered phase including its 
geometrical packing, the free energy functional described here will allow us to study various 
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phenomena of ordered phases. Our work on freezing of simple liquids into crystalline solids 
is in progress and the results will be reported in near future. 
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